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High Energy Laser Beam Propagation in the Atmosphere: 

The Integral Invariants of the Nonlinear Parabolic 
Equation and the Method of Moments 

Robert M. Manning 

National Aeronautics and Space Administration 
Glenn Research Center 
Cleveland, Ohio 44135 

Abstract 

The method of moments is used to define and derive expressions for laser beam deflection and beam 
radius broadening for high-energy propagation through the Earth’s atmosphere. These expressions are 
augmented with the integral invariants of the corresponding nonlinear parabolic equation that describes 
the electric field of high-energy laser beam to propagation to yield universal equations for the 
aforementioned quantities; the beam deflection is a linear function of the propagation distance whereas 
the beam broadening is a quadratic function of distance. The coefficients of these expressions are then 
derived from a thin screen approximation solution of the nonlinear parabolic equation to give 
corresponding analytical expressions for a target located outside the Earth’s atmospheric layer. These 
equations, which are graphically presented for a host of propagation scenarios, as well as the thin screen 
model, are easily amenable to the phase expansions of the wave front for the specification and design of 
adaptive optics algorithms to correct for the inherent phase aberrations. This work finds application in, for 
example, the analysis of beamed energy propulsion for space-based vehicles. 

1.0 Introduction 

When delivering large amounts of power through the Earth’s atmosphere via millimeter or infrared 
beams (i.e., laser beams or beams formed at the output of a millimeter wave antenna system), many 
propagation mechanisms must be addressed that may be potentially deleterious to such power 
transmission. The most obvious one is the ever-present random variation of the atmospheric refractive 
index due to local temperature variations known as ‘turbulence’. This naturally occurring phenomena is 
driven by thermal convection of heat from the Earth’s surface; once the resulting air motion exceeds a 
critical value of velocity, laminar flow essentially evolves into turbulent flow and fluctuations in the 
temperature distribution becomes statistically random (Ref. 1). These temperature fluctuations then act 
directly on the prevailing refractive index, thus rendering the refractive index a random quantity. These 
refractive index variations randomly focus and defocus the intervening electromagnetic wave field. Thus, 
the atmosphere can be considered to be composed of Tenses’ of random focusing and defocusing 
characteristics that, due to the gross atmospheric motion due to wind, move across the beam. This gives 
rise to many beam quality variations; the major ones being beam broadening and beam steering. The 
statistical analysis and modeling of this type of atmospheric propagation as a long and rich history and 
has resulted in analytical descriptions for the impact of turbulence on the operation of systems relying on 
such beam propagation. Many models and descriptions exist for the ‘engineering analysis’ of the 
operation of transmission systems that rely on the propagation of electromagnetic beam propagation in the 
atmosphere (Ref. 2). 

The scenario discussed above may be considered as ‘passive’ electromagnetic wave propagation, i.e., 
the wave field moves through an atmosphere the refractive index of which is determined by other sources, 
not the field itself. However, as the energy density of the beam increases, absoiption of the beam energy 
by atmospheric gas components results in local heating of the atmospheric which does indeed act directly 
on the refractive index causing it to decrease in value . This thermal change of the refractive index field 
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then acts on the electromagnetic wave field causing it to also change, and so on. The propagation scenario 
now becomes an ‘active’ one, whereby the propagating field modifies the very medium it which it exists. 
This heating process is called ‘thermal blooming’ and substantially differs from that of the passive 
propagation discussed earlier. Here, a ‘thermal lens’ is created within the atmosphere by the heating due 
to the energy density of the beam. This ‘self-action’ of the beam will not only bend the beam into regions 
of higher refractive index (beam steering), but convection within the atmospheric fluid will also arise 
which is the source of self induced turbulent flow of the medium. The situation is further complicated 
when one includes the effects of atmospheric wind and aerosols and the abovementioned passive 
propagation effects. Defocusing and other such associated nonlinear thermal blooming distortions of the 
beam cross-section will then result. In extreme cases of very large energy densities, the Kerr effect will 
arise and the propagating beam can subsequently break up into smaller beams, or filaments, which 
severely constrains the amount of energy density that the beam will be able to possess as it travels through 
the atmosphere. Unlike the situation of passive propagation, the thermal blooming mechanism introduces 
nonlinearities into the analysis of the phenomena that substantially complicates a complete mathematical 
description. Complete analyses of these types of propagation scenarios can only be done numerically, 
which was a major activity within the United States and Russia in the late 1980s. Other than the usual 
‘order-of-magnitude’ estimates using the equations of fluid mechanics and wave propagation, only 
numerical modeling of the effects of atmospheric thermal nonlinearities abound in the literature. 
Analytical treatments appropriate for an engineering analysis and assessment of atmospheric propagation 
systems encountering thermal blooming have been lacking and those that do exist hold only for the 
geometrical optics case (Refs. 3 to 5). 

Two of the more important performance parameters of a high-energy laser-beaming scenario are the 
accuracy to aim the beam in the desired direction and the size of the beam at the intended target. In most 
cases, the beam will be focused at the target so the relevant quantity is the size of the focal spot on the 
target. It is the puipose of this memorandum to derive expressions for the beam deflection from the 
intended axis and the broadening of the beam focal spot at the target using only the information that can 
be gleaned for the structure of the nonlinear parabolic equation that describes the propagation process in 
the presence of thermal nonlinearities without actually solving the equation in its entirety. This will be 
accomplished by invoking the integral invariants that exist for the nonlinear parabolic equation as well as 
employing the method of moments to express the statistical moments of the electric field that are used to 
define the deflection as well as the radius of the laser beam. Once this has been done, a thin screen 
approximation solution will be given for the propagating wave field of the laser beam through the Earths 
atmosphere and these results will be applied to derive analytical and graphical results for the beam 
parameters focused to a target above the atmospheric layer. This thin screen model will allow for the 
design and specification for the adaptive correction of the phase aberrations that prevail in such 
propagation situations, and the expressions for the deflection and radius will allow for the evaluation of 
their corrected values. 

2.0 The Moments of the Laser Beam Energy Density and the Nonlinear 
Parabolic Equation 

i i2 

Given the energy density |f?(jc,y,z)|~ of the electromagnetic field distribution of a laser beam in a 

plane r ± = xx + yy transverse to the direction of propagation at a distance z from the origin, one can 
define the mn th moment of the deposition of this quantity in the transverse plane, viz, 


Y , m (z)= J J x' n y"\E(x,y,z)\~dxdy 


( 1 ) 
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Thus, Yoo(z) is the total beam energy at the distance z, Yio(z) is the first moment or centroid of the beam 
energy from the z-axis in the x-direction, y 01 (z) as that in the y-direction, y 20 (z) as the second moment of 
the beam energy in the x-direction, etc. Hence, using these spatial moments of the energy density in a 
plane, one can now form the vector quantity 

4l(z) = — Vt(Yio(z)* + Yoi(z)y) ( 2 ) 

Too (z) 

that can be defined as the normalized overall transverse displacement of the beam energy from the z-axis 
(intended direction of propagation) at the distance z, and 

a eff( Z )= Zt(Y20 ( Z ) + Yo2 (z)) (3) 

Too (zj 

that can be defined as the square radius of the spread of the beam energy about the z-axis. This particular 
moment can be centralized by forming the quantity 

r tf(z) = a 0(z)-K(z) ( 4 ) 

which now gives the radius of the beam energy (i.e., the beam radius) about the axis shifted by A ± (z) 
from the intended direction. 

The analysis of these quantities can be a rather difficult problem in the full nonlinear case of high 
energy laser beam propagation in an absorbing medium which gives rise to a host of thermal 
nonlinearities collectively known as thermal blooming. It thus becomes of interest to ascertain just how 
much of the behavior of these performance parameters can be gleaned from the fact that the equation 
governing the electric field of the laser beam in an absorptive nonlinear medium is given by 

2ik — + ViE + k 2 e T T(\E\ 2 )E = 0 (5) 

dz V ' 

without fully solving specific problem scenarios. Here, V]_ is the Laplacian taken along the direction r 
transverse to the direction of propagation z. Equation (5) is the well-known parabolic equation applied to 
a propagation medium whose variation of permittivity Ae is governed by the change in temperature T of 

the medium from its nominal value which is determined by the energy density \E\ of the field, i.e., 

A £ (r,|E| 2 ) = £ r r(|E| 2 ) (6) 

where £ r = de/dT is the variation of the permittivity with respect to T; in the case of the Earth 
atmosphere, Z T < 0 . The functional dependency of the temperature variation on the field energy density 
is given by the application of the conservation of energy to the atmospheric heat budget, given in its 
entirety by 


P c p 


dT 

lh 


+ V-VT 


= kV~T + a- 


871 1 


( 7 ) 


where p is the density of the atmosphere, C p is its specific heat at constant pressure, K is its thermal 
conductivity, V is the atmospheric wind velocity, and a is the absorption coefficient. (Strictly speaking. 
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one should also include the set of Navier-Stokes equations in the Boussinesq approximation to account 
for the viscous and Archimedean forces that occur during the thermally induced motion of the 
atmosphere. However, for the puiposes of this discussion and subsequent analysis, a description at this 
level is not required.) The implications of the parabolic equation of Equation (5) applied to this nonlinear 
problem will now be investigated for the evolution of the moments of the beam parameters discussed 
above. 


3.0 The Invariants of the Nonlinear Parabolic Equation 

Two integral invariants (sometimes called the ‘integrals of motion’) of the differential equation of 
Equation (5) can easily be derived (Refs. 6 to 8). First, multiplying Equation (5) by E and subtracting 
from it the complex conjugate of the same equation gives 


( 


2 ik\ 




£ *A +e ^: 

dz dz 


+ E*V 2 i E-EV 2 i E* = 0 


( 8 ) 


Using the fact that 



EE" and simplifying where possible yields 


d E ^ 

2 ik- + V ± -(£*V 1 £-£V 1 £*) = 0 (9) 

dz v ' 

Finally, integrating Equation (9) in the transverse plane yields and applying Gauss’s Law to the second 
member of the relationship and letting the surface approach infinity, at which the surface terms approach 
zero, one finally obtains 



E\ dxdy = const = y 00 


( 10 ) 


which is just a statement of the conservation of energy for any transverse plane. (Absorption of the wave 
field is present which, after all, accounts for the heating and subsequent nonlinear action of the medium. 
However, the action of the absorption on the overall energy density of the laser beam itself has been 
neglected in Equation (5). Incorporating it into the subsequent development is straightforward but will not 
be done here.) 

Second, multiplying Equation (5) by dE"/dz and adding to it the complex conjugate of the same 
equation gives 


^ Vi E+ d ^- Vi E* + kh T T (|£| 2 )^ = 0 (11) 

As done above, one now integrates this relationship over the transverse plane. Each of the first two terms 
on the left side of Equation (1 1) can then be simplified by first using Green’s first identity and letting the 
surface terms vanish for the same reason as employed above: 

OO OO 0l^r|2 

dxdy + k 2 e T J J r^|E|~ j ^ d xdy = 0 (12) 

Defining the quantity G^|Z?i j such that 


-j j 9 X ~ V 
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and substituting it into Equation (12) and simplifying yields 


d_ 

dz 


j J' v ' v ' 


\dxdy + k 2 e T — 
' ' dz 


J J G^\E\*\dxdy = 0 


(14) 


which becomes 



V,£-V,£* 


-k e t G 


(id 3 ( 


dx dy 


= const = H 


(15) 


which is essentially the Hamiltonian for the propagating wave field. 


4.0 Conservation Laws for the Nonlinear Parabolic Equation 

The two integral invariants given above should now be augmented with the conservation laws that 
prevail quantities that issue from Equation (5). A streamlined method to accomplish this is to employ 
field theory and find the Lagrangian corresponding to Equation (5) and from it obtain the energy 
momentum tensor from which the conservation quantities and laws would follow (Ref. 9). However, a 
direct approach will be used using Equation (5) directly. To this end, the first such conservation law was 
already derived and is given by Equation (9). Defining 


^=-^(£*v ± £-EV ± zr) 


which is the energy flux vector of the field, one has 


9|E| 2 

dz 


+ V p -5 X =0 


(16) 


(17) 


This is one of two prevailing conservation laws. 

Consider now the other equation which governs the energy flux vector 5 ± ; the idea here is to obtain 

— . .2 

an equation for S ± analogous to that of Equation (17) for \E\ . Using Equation (16), 


f i ] 

dE* - 

' dE ' 

dE - * - 

( dE * T 


V,£ + £V 


V,E -EV, 


,2k, 

dz 1 1 

l dz J 

dz 1 1 

l dz )_ 


Employing Equation (5) and simplifying gives 


^ = f— 1 [ -K*F + E*W,K-KF*+EW,K*] 
dz {2k ) L 1 1 J 


(18) 


( 19 ) 


where 
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( 20 ) 


F = V ± E, K = V ± -F + k 2 e T TE 


Making use of the identity 

w ± (e*k) = f*k+e*w ± k 

as well as its complex conjugate allows Equation (19) to become 

f « 


( 21 ) 


95 x 

dz 


— ) \-2FK* -2F*K + V 1 (e*k) + V ± (eK*) 

2 k J 


Two more identities are now needed: 


V; \E\ 2 =2F-F + E*y , -F + EV : ■ F* 


and 


( FW l F* + F*V ± f) =—(f i F- ) + —( E*F ) - — (e ■ E * ) 

1 1 ± fi 9x/ ' 1 ’ dxj ' J > dxE ’ 

(The latter comes from the fact that F x V x x F = 0 since F = V X E , etc.) After a bit of algebraic 
manipulation, one finds that 


( 22 ) 


(23) 


(24) 


'95, ^ fO 2 


v dz J 


\2kj 


- 2 -(^' + ^) + -( v i|£| 2 ) 


+2fc 2 £ T 


f 9 / | ,2\ 

9 . a] 


— \T\E\\ 

-T \E\ 


9x V 1 1 / 
( ax j 

?)x j J 

i j 


(25) 


Finally, using Equation (13) and simplifying, Equation (25) gives 

f o? 5 7 i 

\k) 


95 x 
v dz j 


dx j L 


(V X E* ), (V ± E) . + 1 (r |E| 2 - G (|E| 2 )) 8, y 


if* 


(26) 


Equation (26) can be written as the tensor relation 

95 


dz 


- V x <2 = 0 


(27) 


where Q is the stress tensor 
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( 28 ) 



r i \ 2 ' 

JJ [4 


i(vi|£| 2 )5, ; -i(v ± £),(v ±£ *). 

(v,£-) ( (v l£ ). + l -e 


l r 


— k e T T \E\ -G 


(l £ l 2 ))5 9 


The conservation laws of Equation (17) and (27), as well as Equations (10) and (15), which prevail 
for Equation (5) are now sufficient to obtain expressions for the beam parameters specified by Equations 
(2) to (4). 


5.0 The Displacement and Broadening of a Laser Beam Experiencing 
Thermal Blooming 

The displacement from bore sight of a laser beam which induces thermal nonlinearities during 
propagation can now be easily calculated using Equation (2) involving the first order spatial moment of 
the energy density. To this end, differentiating Equation (2) (noting the constancy of y 00 by Equation 
(10)), and using Equation (17) yields 


d ^ L = J J Hv x -5 ± )x + y(V ± -S^dxdy 

Integrating by parts and noting that the surface terms approach zero at infinity gives 

^ = ^ii S±dxdysAiz) 

This is sufficient to describe the displacement in that 

c lAdA = J_ j j ^dxdy = — f f V ± • Qdxdy = 0 
dz 2 Y(K) dz Y 002 L 2 L 


(29) 


(30) 


(31) 


upon using Equation (27); the last term vanishes upon using Gauss’s Theorem and, once again, letting the 
fields and their derivatives vanish at infinity. Thus, A(z) = A = const and the overall beam displacement 
is given simply by integrating Equation (30), i.e., 

A ± (z) = Az (32) 

where it is assumed from hear on that A x (0) = 0 , i.e., the beam is initially pointed down its intended 
axis. Thus, the beam deflection from its intended axis is a linear function of distance travelled. 

Similarly, the associated overall beam broadening can be found by first differentiating 
Equation (3) and using Equation (17), 

— ~j~ = — ~ J °\(x 2+ y 2 )(V±-Si_)dxdy (33) 


Integrating by parts gives 


NAS A/TM— 20 12-217634 


7 



da fz) 2 r r - / .. 

— = r ± •5_ l dxdy = B(z) (34) 

dz Too jLjL 

Differentiating once again and using Equation (27) yields 

eff( z ) = — 1 j r ± -(V_l Q)dxdy (35) 

dz Too _ M -oo 

Integrating by parts and letting the surface term vanish, 

] ff 2 ^ = — — J °\[h-Q)dxdy = — — J j ^ dx dy (36) 

dz Too 2 L -00 Too tL-L i=i 

where « is the surface normal. In order to assure that this result is a constant and that no further 
derivatives are needed, one must evaluate the required trace of the tensor Q . Using Equation (28), one 
has 



1. 

2 


\k) 


2 

| V_ l -(£*V_ l E + £V ± E*)-2V_ l E 


■W L E* + 2k 2 e T 




(37) 


Introducing this into Equation (36) and using the Gauss Theorem to rid of the divergence term finally 
gives 


d 2 a 2 e ff (z) 

dz 2 


/j_Y 
vToo A 


JJ[ 


-2V , E ■ V , E* + 2k 2 e T T \E\ - 2 k 2 e T G 


(id 2 )' 


dx dy 


(38) 


This expression can be stated in terms of the invariant of Equation (15), viz, 


d 2 a 2 ff (z) ( i V,V 


dz" 


1_ 
(.Too 


]_ 

\kj 


2H —2k 2 £■, 


j Jt 7 ' id 


2 -2g(|eY 


dx dy 


(39) 


Hence, this expression is constant so long as the integrals vanish within the brackets. Thus, in the case of 
the Kerr effect, from Equation (6), one has that the temperature is replaced by T => A\E\ where A is 

some parameter independent of the electric field intensity, then by Equation (13), G (|Yj = ( 1/2 ) A | /v| 4 

and the integrand is identically equal to zero. However, in the case of thermally induced nonlinearities, 
the temperature field is related to the electric field intensity through Equation (7). In the case of forced 

convection due to wind and where large Peclet numbers prevail (see Appendix A), dT /dx = P|/r (x,y,z)| - 

(for a wind velocity in the x-direction), and this will give rise to much more complicated functional 
dependencies. A general analysis of this problem is beyond the scope of this particular work but as shown 
in Appendix A, in the case of forced convection, 


jj T|£| 2 -2 g(|£| 2 ) 


dx dy~0 


(40) 
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Hence, 


d 2 a 2 eff (z)_ 

( 1 7 


dz 2 

vT00 y 

\k) 


2 H = const = C 


(41) 


Integrating this expression and remembering Equation (34) gives for the overall effective beam radius 


1 




(z) = E Cz 2 + B(0) z + a 2 (0) 


(42) 


and from Equations (4) and (32), for the effective radius about the deflected axis, 

7 \z) = Uz 2 +B(0)z + a 2 ^(0)-(Az) 2 


' eff ' 


(43) 


Hence, Equations (32) and (43) give the universal behavior of the deflection and broadening of a laser 
beam governed by Equation (5) without solving the entire propagation problem. The constants A , B(0), 
and C are related to the beam parameters via Equations (30), (34), and (41), respectively. 

As a quick check on the validity of Equation (43), consider a Gaussian beam wave propagating in free 

space. Thus, T ^|£ , |" j = 0 in Equation (5) and a solution of the resulting equation can be given in the form 


E(f,z) = - 


A) 


-exp 


ka d A 


1 + ia d z 


1 + ia d z 

where A 0 is the initial maximum amplitude of the wave in the output aperture and 

2 


(44) 


OCj — OCi “I - iOCo ? OCi — 


kW,t 


F 


(45) 


with W 0 being the waist size of the beam as it leaves the aperture (the effective radius is a eff (0) = W 0 /y[2 

) and F being the focal length of the beam. (The beam parameters a,,, ex, , and a 2 ai 'c not to be confused 
with the absoiption coefficient a of Equation (7).) Hence, one has 


where 


(46) 


1 p|2 4 2 

k<x l 2 

\ E \ = . f yexp 

A(z) 

[ AW" J 


A(z) = (l-a 2 z) / +a 2 z 2 


and from Equation (10), 


Too 


(47) 


(48) 


which is the total energy contained within the beam. Employing Equation (16) gives 
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{xx + yy) 


(49) 



V 


(oq + a 2 )j 

A(z) 


Equation (30) then yields 


and from Equations (34) and (41) 


A = 0 


(50) 


f 2 A 


V ka i J 


B(z) = ~ — (a 2 -(ar + a;)z 


C = 


2[a\ + a 2 ) 


ka , 


(51) 


(52) 


Using these values in Equation (43) finally gives for the effective radius of the beam in a transverse plane 
at a distance z from the source aperture 


4(4='f[( | -a 2 4+<* 2 ] = ^A(z) 


r ' A l <• ( 53 ) 

which is the value found using from Equation (46) by equating the exponential argument to - 1 . 


6.0 Application of the Foregoing to High Energy Laser Beam Propagation 
Through the Atmosphere 

The displacement and widening of a very high power laser beam due to atmospheric thermal 
nonlinearities are two very important performance parameters in the case of transmission of, e.g., laser 
energy from the Earth’s surface to power a spacecraft in a low Earth orbit (LEO) for subsequent 
propulsion into, e.g., a geosynchronous orbit (GEO). In this instance, a thin screen approximation can be 
made whereby only a very small length L of the entire longitudinal propagation path distance D is 
intersected by the atmosphere, i.e., L « D. To simplify matters, consider instead of the field E, its 
amplitude and phase representation 


E(r,z) = e(r,z)exp[-i§(r,z)-\ (54) 

Substituting this into Equation (5) and separating the real and imaginary portions of the resulting equation 
gives 


2k ^ ^ + (' V J)) 2 + k 2 e T T (|E| 2 ) = 0 
oz e v ' 1 ' 

2k - — 2 V j_c() - V ± c-(v5_{|))c = 0 
dz x ' 

Similarly, using Equation (54) in Equation (16) gives 


(55) 

(56) 
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(57) 


Sj_=~ 


e 2 V§ 


and using this intermediate result in Equations (30), (34), and (41), one has 

1 (■ 


A = — - 


Too 


\ J J k| 2 (V ± (l ))dxdy 


(58) 


B = — 


Too 


i) J J \ e \ 2 ( r ±-v±$) dxd y 


(59) 


C = 


f 2 V i 5 2 


woo 


r 'JJ IM 2 + |c| 2 ( VJ)) 2 - k 2 e T G(\E \ 2 ) j dxdy 


where from Equation (10), 


(60) 


Too 


= J J \e\ 2 dxdy 


(61) 


The system of Equations (55) and (56) can be considerably simplified by using the fact that the layer 
through which the laser beam will encounter thermal nonlinearities is very thin with respect to the entire 
propagation path. Hence, letting only the phase of the field be perturbed within this thin layer and, thus, 
neglecting amplitude effects, one must solve the phase equation across this thin atmospheric layer, i.e., 


+ (v x ^) 2 + £ 2 £ x r (|£| 2 ) = 0 


(62) 


where the perturbing temperature field T ,z)|~ j is given by Equation (A.3) of Appendix A. As 
shown in Appendix B , one has in this “thin screen” approximation 

4> = — - — — y^(a 2 ~{a\ + al ) z) + VrcO (P 0 ) exp (- tojy 2 ) ( 1 + erf ( 

2A(zj 


(63) 


where 




y2 kW 0 j 


>0 


(64) 


which convolves the initial power P 0 within the laser beam and the power P T that is characteristic of 
atmospheric thermal effects, viz, 


P T - ,2 


\[2npVC 


k W 0 e T a 


( 65 ) 


The negative sign occurs in the definition of Equation (64) since, as noted earlier, e r < 0. 
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The evaluation of Equations (58)-(60) commences with Equation (61); by using the first relation of 
Equation (B.3) of Appendix B evaluated at the output plane of the laser in Equation (61) and performing 
the indicated integrals in plane polar coordinates gives the same result as Equation (48), i.e., 


Too = 


ka l 


( 66 ) 


The calculation of the vector A requires the evaluation of the transverse gradient Vj_(j) . Using Equation 
(63) one has 





{xx + yy) 



exp 


(-toi/) 


(67) 


Substituting this expression into Equation (58) as well as that of Equation (B.3) evaluated at z = 0 (since 
A = const by Equations (30) and (31) and so can be evaluated at any longitudinal location), one finds that 
all terms involving first powers of x or y vanish upon integration in the transverse plane leaving 



Hence, by Equation (32), the beam is deflected by 



indicating that the beam moves against the atmospheric wind direction. This phenomenon is easily 
explained by the fact that the beam wants to move into regions of higher refractive index, i.e., into cooler 
regions about the beam axis. 

Since the beam deflection adds an asymmetry of the beam along the x-axis, it proves to be beneficial 
to separate out the x and y contributions of the remaining B and C parameters of Equations (59) and (60). 
To this end, Equation (59) can be written 


B = B+ B, 


where 


B= — 


Too 


- j J J \e\ 2 (xx-V^)dxdy 


B=~- 


Yoo 


“1 J J \e\ 2 (yy-V ±§)dxdy 


(70) 


(71) 


( 72 ) 


Using Equation (67) and once again noting that terms containing first powers of x and y vanish upon 
integration, one obtains 
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One must now separate Equation (60) into x and y component contributions. Thus, without loss of 
generality, one can write 


c = c x +c y 


(74) 


where 


r = 


( 2 V^ 2 


l Too 


— j j J J |(x-V ± c)| +\e\ |x-V ± (|)j dxdy — Q J J k 2 z T G{^E\ j dxdyl (75) 


and 


c,= 


' 2 Vl' 2 

v Too 


J J |(t • v _l^)| + k| 2 (v - V_l4>) _ dxdy-(l-Q) J J ^ 2 8 r G||E| 2 ) dxdy 1 (76) 


where the j terms have been separated out for convenience. It is very important to note that 

G ^|£|- j cannot be arbitrarily apportioned between the arbitrary x and y components of the coefficient C. 

The j term has its origins with the conservation of energy given by Equation (15) and is thus 

associated with the deposition of energy within the propagating laser beam. The beam deflection adds 
asymmetry to the beam and energy will be deposited along the x and y directions differently. In what is to 
follow, a fraction Q of the energy density will be taken to correspond to that along the x axis and the 
remaining 1 - Q fraction along the y axis. An approximation procedure will be given below to determine 
the division of energy among these components and to find a value for Q. 

Consider, at the outset, the evaluation of the first integrand terms of Equations (75) and (76). Once 
again, working in the z = 0 plane, one has from the first of Equation (B .3) 


(x • V ± <?) = A^k 2 a 2 x 2 exp (-/coq (x 2 + y 2 ) j 
and a similar relation for the y component. One thus has 

J J|(* = J dxdy = 


(77) 


(78) 


Using Equation (67), one also has 

. \2 


e\ (x-V(|)) = A,J exp^-£oq (x 2 + y 2 j ( ka 2 x )“ + 4/ca 2 xO ( P (l ) 


•exp 


[-ka^x 2 j exp (-&oc,y 2 j + 40 2 (P 0 ) ka x • • exp (-2fcoqx 2 j exp (-2£oqy 2 j 


( 79 ) 


and similarly, 
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exp^-koq (x‘ + y 2 )jy 2 ( ka 2 ) -4k 1 a 1 a 2 \fn^>[P 0 ) 

■ exp [~ka x y 2 j ^ 1 + erf^koqxjj + +4k 2 cx 2 7t0 2 (P 0 ) exp [~2ka x y 2 1 + erf^koqxjj 


These yield 


J Jh 2 ( 3 


dxdy = At j - + Aq 2 | TtO 2 (P 0 ) 


J J M 2 (f Vc^) 2 ^y = 4^-A 0 2 7i^^O(P 0 ) + A 0 V /2 -^O 2 (P 0 ) 


The evaluation of the second integrand term in Equation (76) commences with the use of 
Equation (40), i.e., 


k 2 e r J J g||E| jdxdy = — k 2 e r j Jr|Zs| dxdy 


One now employs Equations (46) and (A.3) in the right side of this relation and obtains in the z - 0 plane 


J J T\E\ dxdy\_ =Q = — (3 J J exp [~ka x x 2 ) exp [-2ka x y 2 ) |l + erf ^ka x x j j dxdy 


2 \ ka i LL 

2 ^ k(x^ j V2 


Using this in Equation (83) and rearranging terms gives 


k\ T J J G[\E\ 2 )dxdy = ^^^ 


It may be noted that Equation (85) is independent of the thickness L of the absorbing medium. This is due 
to the level of estimation that the ‘thin screen’ approximation yields; in that approximation, only the 
phase and not the amplitude of the field intensity of the beam is perturbed. It would only be with the next 
approximation beyond that used in Appendix B that would add a perturbation to the amplitude of the 
field. This limitation must be kept in mind in the applications. 

As discussed above, the division of this quantity among the x and y components cannot be arbitrarily 
done. The value of Q, introduced above, will be calculated shortly. Thus, substituting Equations (66), 
(78), (81) and Equation (85) into Equation (75) results in 


c '=t[ i+ | + ^ 2 < p 4J+^J*t.) 


2 \kL 


and, similarly, putting Equations (66), (78), (82), and Equation (85) into Equation (76) gives 
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(87) 


r _^L 

• v k 


( 


+(1 -Q)Ji ir H p o) 


a 1 


V 


ar 


a, 


2 TC 


3^3 


kL 


In applying Equation (43) to calculate the beam radius about the local axis, care must be taken in the 
specifications for the B and C parameters using Equations (73), (86) and (87). In fact, the derivation of 
Equations (73), (86) and (87) require the corresponding parameters in Equation (43) to be multiplied by a 
factor of 2, i.e., 


^,Az) = ^2C x )z 2 +{2B x (0))z + a^(0)-(A.xzf (88) 

and 

4 rW = ^( 2C > 2 + ( 2 fi r ( 0 ))^ + 4(°)-( i -^) 2 ( 89 ) 

Using Equations (88) and (89), one has from the foregoing, after rearrangement of factors and 
simplification, 


a i^ r eff,x i z ) - 


'eff,x 

Wr 


(l-a^z) 2 +a 2 l + -O 2 (P 0 ) z 2 +Qa l 
V 3 



m 


and 


«1 kr eff,y (Z) = - Jf = ( l ~ + 


W f 


1 + 


16 n 

~9~V3 


® 2 { p o) 


f 1 \ 


+(l-G)a 1 ^^-Jo(P 0 )r +2a 1 ^O(P 0 )z(l-a 2 z) 


(91) 


The case of a focused beam wave, i.e., F = D will be considered, without addressing the methodology 
with which to achieve it in practice. In this instance, the above relations reduce to 


% =a4i4®qp 0 )W + eaJ^U®( p o)° 2 

Wo V 3 4 V 2 ( L ) 

foe 


f 2^ 


(92) 


and 


'eff,y 

W n 7 


= ar 




foe 


( 1 \ 


(93) 


An expression for the value of Q can now be obtained. The first member of the right side of Equation (92) 
is the semi-axis of the beam in the x direction and that of Equation (93) is the semi-axis of the beam in the 
y direction. As a first approximation, it can be assumed that the ratio of the lengths of these semi-axes is 
equal to the corresponding ratio of the energy deposited along these axes, i.e., 
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l+3° 2 ( P o) 


Q 


l + ^^O 2 ^) 1 Q 

9 S 


( 94 ) 


Solving this expression for Q yields 


La, 


Pt+KJo , 

Q = Q{P 0 )= 4 


2P;- + {K l+ K 2 )P^ 


(95) 


where K 1 = 5/3 and K 2 = I (m/9xf?> . This expression is well behaved at the limits lim Q = 1/2 

(indicating that the energy is deposited equally along the axes when no thermal action is prevalent) and 
lim Q = 0.341 . Once again, it must be remembered that this approximation procedure is in lieu of a 

P 0 ^°o 

second-order solution to the equations of propagation in which the field amplitude of the beam wave is 
also given as a function of propagation conditions along the beam. 

Numerical application of the foregoing to focused beam propagation in a nonlinear thermal medium 
commences with the rewriting of Equations (92) and (93) in a dimensionless form. Additionally, the 
expression for the corresponding beam deflection A x as given by Equation (69) will be treated. First, 
dealing with the latter expression, Equation (69) and be rewritten in the dimensionless form (deflection 
from intended axis normalized by the initial beam radius) 


A„ = 


W n 


= V2JLo(P 0 ) * /,O(P 0 ) 


kW r 


V2 


(96) 


where l dD = 2 L)/ kW (2 is the diffraction length associated with the aperture size and the total propagation 
distance D to the target. Figures 1 and 2 depict the variation of the normalized deflection of the beam 
versus the beam power in kilowatts, for l dD = 0.02 and l dD = 0.2 , respectively; these correspond to 
realistic beam and propagation conditions. The direction of the deflection is along the -x axis, i.e., it is 
directed into the wind where regions of higher refractive index prevail. As heat is removed from the 
propagation channel along the x axis, it continues to accumulate along the y axis of the channel. This 
gives rise to an asymmetry in the beam cross-section as is shown by considering the beam radii. 

Similarly, for the beam radii at the intended focus, one can write Equation (92) as 


4..M? 


(a,D) 


= r F,x = 


l + -^o) 


A 


foe 


+ e,if 


( 


a X L 


®(?o) 


(97) 


The left side is essentially the square of the ratio of the beam radius in the v direction normalized to the 
intended radius of the focal spot on the target, 2 D/kW 0 . The same for Equation (93), i.e., 


r F,y = 


16 K o / \ / \ 71 

Ti/f ( °)J + ( 1_e )V2 


ajL 




(98) 
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In both of these relations, the quantity IdL = a l L is the diffraction length with respect to the aperture size 
and the thickness of the atmospheric layer L. Figures 3 and 4 depict the beam radii variation versus beam 
power for two representative values of IdL ■ 



Figure 1 . — Beam deflection from intended axis normalized to initial 
beam radius versus beam power in kilowatts for l dD = 0.02. 



Figure 2. — Beam deflection from intended axis normalized to initial 
beam radius versus beam power in kilowatts for l dD = 0.2. 
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Figure 3. — Beam radii along x and y axes normalized to focal spot 
radius at target versus beam power in kilowatts for Ijd = 0.001. 



P 0 , KW 

Figure 4. — Beam radii along x and y axes normalized to focal spot 
radius at target versus beam power in kilowatts for I ( id = 0.01. 


In addition to the evaluation of the deleterious effects of the atmosphere on the propagation of high- 
energy laser beams, the methodology introduced above is amenable to the analysis of the design and 
specification of adaptive optics approaches whereby the phase front of the transmitted laser beam is 
appropriately modified to mitigate the atmospheric effects (Ref. 10). For example, due to the elliptical 
profile of the radiation in the transverse plane, the use of the well-known Zemike polynomials (Ref. 11), 
which are defined on a unit circle, are not appropriate for this application. Instead, the transmitted phase 
front <\> T (x,y) at the output aperture of the transmitter can be expanded in a set of orthogonal 
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polynomials associated with a Gaussian weight function, i.e., the Hermite polynomials H„(x) and H m (y). 
Thus, 


M N 

4>r = Y,Yj a n>n H n (4 H m (?) (") 

n m 

The expansion coefficients a nm are determined by applying one of several ‘performance metrics’. For 
example, for power beaming applications, it is desired to shape the phase front so as to minimize the 
radius of the beam at the target. Such a representation is straightforward to use in Equation (63) and 
subsequently, in Equations (71), (72), (75) and (76) to obtain the corrected radius values for the beam. 

7.0 Summary and Conclusion 

The structure of the nonlinear parabolic equation that describes the propagation of a high energy laser 
beam through a medium with thermal nonlinearities is exploited to describe the moments of the 
associated electric field in order to capture the universal behavior of the displacement and radius of the 
beam; these are two of the more significant performance parameters for the operation of the high energy 
laser transmission system. Once this was done, a simple thin phase screen propagation model of the Earth 
atmosphere with a wind was advanced and applied to the equations for the beam displacement and radius. 
The presence of the wind of course broke the symmetry of the beam into one with an elliptical cross 
section. Expressions were then developed for the behavior of the semi-axes of the beam for the 
propagation path out of the atmosphere. Graphical results were then presented for the deleterious action of 
the atmosphere on the propagation of the beam. Since only phase perturbations were admitted into the 
model, all the effects placed on the propagating beam are phase based and can thus be mitigated at the 
transmitter. The elliptical cross-section of the beam dictates that the circular Zemike polynomials are not 
sufficient to describe the phase front and that a Hermite polynomial based phase decomposition be used. 
Such a model was exploited for the analysis of laser power beaming for space propulsion (Ref. 12). 
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Appendix A 

Equation (7) can be simplified straight away by first specializing the wind velocity to be along the x 
axis of the coordinate system, V = Vx , where x is the unit vector. Further, the diffusivity V 2 T will be 
taken to have contributions only in the direction transverse to the laser beam propagation, i.e., 

V 2 r ~ V 2 T . Finally, only the stationary heating case will be considered for this particular - discussion 
whereby the time derivative can be dropped. (That is, the CW radiation is taken to be acting long enough 
for the steady state case to be achieved.) Applying these considerations to Equation (7) allows it to be 
written 



XVlT = 


i £ |2 

87tpC p 


(A.l) 


where % = K/p C p is the associated thermal diffusivity. There are two heat transfer mechanisms that are 

described by this relation, i.e., forced convection due to the wind given by the first term on the left side of 
the equation and diffusion by the second term on the left side. The order of magnitude of the 
corresponding derivatives are dT /dx ~ T '/r e g and V\T ~T j Q j where r eff is the effective laser beam 

radius. Consider now the ratio of the coefficients of the two terms within the brackets of Equation (A.l), 
i.e., t) = Vr e ff/x . This dimensionless ratio is known in fluid mechanics as the Peclet number. Taking the 

nominal wind velocity V = 4.4 m/s (10 mi/hr), r eff = 10 m and using the documented value for the thermal 
diffusivity of the atmosphere % = 2.12 x 10 -5 nr/s, one has that T| = 2.1 x 10 6 » 1 allowing one to 
neglect the second term in the brackets. In this instance, Equation (A.l) becomes 

^ = (3|E| 2 (A .2) 

dx 


where (3 = ac j {%iiV pC . This equation can now be easily integrated using Equation (46). The result is 


X 

r||£'(r,z)| 2 j = (3 J \e(x' , y,z)\~dx' 


= P 


K 


n 


\ 1/2 


kaiA(z) 


exp 


tai o 

;y 


A(z)‘ 


( 


1 + erf 


V 


fka^ m 

v A ( z )y 


(A.3) 


jy 


Now in order to substantiate Equation (40) for this case, one proceeds to define 


/(r(|E(x)| 2 ),G(|E(x)| 2 ))= J j[T\E\ 2 -2G(\E\ 2 )y x dy 


(A .4) 


and attempts to show that the expression is equal to zero. To this end, one forms the derivative 


dl 

dx 



d\E\- dG 

+ - ! — L r-2 

dx dx 


dxdy 


(A .5) 
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which, if also zero, can be used to argue that Equation (A.4) vanishes do to conditions at x — > °° . From 
Equation (13), 

3 M 

G( Z{x)) = J T(Z')dZ', Z(x) = \E{xf (A.6) 

o 


Differentiating Equation (A.6), 


^ = ^hir(=) + f^i/=' 

dx dx ' dx 


(A .7) 


substituting into Equation (A.5), using Equation (A.2), performing the indicated functional integration, 
and simplifying finally gives 


dl 

dx 


=-JJ 


d\E - 

dx 


dxdy 


(A.8) 


which does not appeal - to vanish. To be sure, proceeding to integrate Equation (A.8) by parts yields 


— = f f |£| — dxdy = -(3 f f| E\ 4 dxdy 

dx J J dx J J 


(A .9) 


I |4 

upon remembering Equation (A.2). The \E\ dependency is reminiscent of the Kerr effect that appears at 

field intensities much higher than for thermal blooming effect. To be sure, in order to secure numerical 
estimates, one uses Equation (46) in Equation (A.9) to obtain 

^ = K^o 2 (A. 10) 

dx 4 


To rid of the amplitude factors and put them into a more familiar quantity, i.e., the total power P 0 within 
the laser beam, one uses the relation 





(n w<n 

l 871 J 

l 2 J 


(A.ll) 


as well as remember the definition of (3 to write Equation (A. 10) as 

df _ 1 6aP ( ; 

dx ~ cW(jVpC p 


(A. 12) 


Using the values of the parameters given above and additionally taking W 0 ~ 1 , this gives 

dl /dx ~ 4 x 10 -5 for P 0 = 1 x 10 6 watts. Hence, dl/dx effectively vanishes for nominal laser power and 
thus Equation (A.4) can also be taken to vanish, i.e., 
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(A. 13) 


jj[^-2 G (| E f) 


dx dy~0 


A more complete treatment would involve a more general solution to the temperature diffusion equation, 
Equation (A.l), but will not be dealt with here. 
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Appendix B 

Writing the total phase as the sum of the unperturbed phase (|) 0 (i.e., defined by Equation (62) with 
T (|£f I = 0 ) and the perturbed component (j) , one has, upon introducing (|) = (J) 0 + tj) into Equation (62), 
the system of equations 

2^ + (V ± 4) 0 ) 2 =0 (B.l) 


2k y + 2 Vj_c)) 0 • V ± ^ + (V x {^) 2 + kh T T (|£| 2 ) = 0 
Now, equating Equations (44) and (54), one finds that 


e(r,z)= 


A 

1 + ia d z 
kr 2 

= 2A {zj 



f ko. x 

exp 

{ 2A(z) 

(«2 

-(a 2 + a: 


(B.2) 


(B.3) 


the second relation of which, indeed, can be shown to satisfy Equation (B.l). 

A solution for Equation (B .2) commences with multiplying it by 1 Ik and considering k — > °° 
(relative to all the other inverse lengths in the problem) which gives 


2^ + fe r r(|£| 2 ) = 0 

dz ' ' 


(B .4) 


One needs to use Equation (A.3) to integrate this expression. To this end, one applies the following ‘thin 
screen’ approximation; For the case considered here, the perturbing atmospheric layer of thickness L, in 
which the thermal nonlinearities take place, is such that L«D where D is the total distance to the target 
between LEO and GEO. Thus, within this region, relative to the total propagation distances involved, one 
can let A(z) ~ 1 for the longitudinal distances z occurring within the atmospheric layer. Thus, one can 
write 


T\\E\ 


j ~ (3 A^ Vk ^ j = exp ^-fcoqy 2 J ^ ^l + erfj^&oqxjj 


A 

2y 


(B.5) 


Substituting this result into Equation (B.4) and integrating across the thickness L of the atmosphere 
yields, upon rearranging factors, 


$ = - 


f n \ 


\ P T J 


\[nL 

j 


( n 


v^y 


exp 


(-fcoqy 2 ^! + erf^ka l x^j 


(B -6) 


where 




'cAf 

r *w 0 a 

v 87I y 

l 2 J 


(B .7) 
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is the initial power of the beam wave of waist radius W 0 and 


V27tpVC„ 


p = : >— 

T k 2 W 0 \e T \a 


is the characteristic power for the onset of thermal blooming. Defining 


V 07 P T 2kWk 


(since e T < 0 ), Equations (B.3) and (B.6) yield for the total phase of the laser beam leaving the thin 
atmospheric layer with a wind of velocity V 

J^l 2 + y 2 i 

<|) = — - — — — ^cx 2 -(a 2 +oc 2 )zj + VttO(P 0 )exp^-to 1 y 2 j^l + erf^tojxjj ( 


(B.10) 
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